Exportemos a R nuestro conjunto de datos para poder comenzar el análisis:
# Exportamos los datos con la función import de la
# librería rio.
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2
datos_sucios<- read_excel("AirQualityUCI.xlsx")
Para mejorar la claridad y evitar confusiones, vamos a proceder a renombrar todas las variables.
# Renombramos las variables para que sea más fácil poder
# trabajar con ellas.
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
datos_sucios<- rename(datos_sucios,
Fecha = "Date",
Hora = "Time",
Monoxido_carbono = "CO(GT)",
Oxido_estaño = "PT08.S1(CO)",
Hidrocarburos_no_metanicos = "NMHC(GT)",
Benceno = "C6H6(GT)",
Titanio = "PT08.S2(NMHC)",
Oxidos_nitrogeno = "NOx(GT)",
Oxido_tungsteno = "PT08.S3(NOx)",
Dioxido_nitrogeno = "NO2(GT)",
Oxido_tungsteno_NO2 = "PT08.S4(NO2)",
Oxido_indio = "PT08.S5(O3)",
Temperatura = "T",
Humedad_relativa = "RH",
Humedad_absoluta = "AH")
Comprobemos si el cambio de nombres se realizó de forma adecuada:
# Veamos si se realizó el cambio en los nombres:
names(datos_sucios)
## [1] "Fecha" "Hora"
## [3] "Monoxido_carbono" "Oxido_estaño"
## [5] "Hidrocarburos_no_metanicos" "Benceno"
## [7] "Titanio" "Oxidos_nitrogeno"
## [9] "Oxido_tungsteno" "Dioxido_nitrogeno"
## [11] "Oxido_tungsteno_NO2" "Oxido_indio"
## [13] "Temperatura" "Humedad_relativa"
## [15] "Humedad_absoluta"
Con las variables renombradas y una nomenclatura más clara, veamos un resumen rápido de los datos:
# Resumen rapido, conociendo los tipos de datos
skimr::skim(datos_sucios)
| Name | datos_sucios |
| Number of rows | 9357 |
| Number of columns | 15 |
| _______________________ | |
| Column type frequency: | |
| numeric | 13 |
| POSIXct | 2 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Monoxido_carbono | 0 | 1 | -34.21 | 77.66 | -200 | 0.60 | 1.50 | 2.60 | 11.90 | ▂▁▁▁▇ |
| Oxido_estaño | 0 | 1 | 1048.87 | 329.82 | -200 | 921.00 | 1052.50 | 1221.25 | 2039.75 | ▁▁▇▅▁ |
| Hidrocarburos_no_metanicos | 0 | 1 | -159.09 | 139.79 | -200 | -200.00 | -200.00 | -200.00 | 1189.00 | ▇▁▁▁▁ |
| Benceno | 0 | 1 | 1.87 | 41.38 | -200 | 4.00 | 7.89 | 13.64 | 63.74 | ▁▁▁▇▅ |
| Titanio | 0 | 1 | 894.48 | 342.32 | -200 | 711.00 | 894.50 | 1104.75 | 2214.00 | ▁▅▇▂▁ |
| Oxidos_nitrogeno | 0 | 1 | 168.60 | 257.42 | -200 | 50.00 | 141.00 | 284.20 | 1479.00 | ▇▆▂▁▁ |
| Oxido_tungsteno | 0 | 1 | 794.87 | 321.98 | -200 | 637.00 | 794.25 | 960.25 | 2682.75 | ▁▇▃▁▁ |
| Dioxido_nitrogeno | 0 | 1 | 58.14 | 126.93 | -200 | 53.00 | 96.00 | 133.00 | 339.70 | ▃▁▇▅▁ |
| Oxido_tungsteno_NO2 | 0 | 1 | 1391.36 | 467.19 | -200 | 1184.75 | 1445.50 | 1662.00 | 2775.00 | ▁▂▇▅▁ |
| Oxido_indio | 0 | 1 | 974.95 | 456.92 | -200 | 699.75 | 942.00 | 1255.25 | 2522.75 | ▁▇▇▃▁ |
| Temperatura | 0 | 1 | 9.78 | 43.20 | -200 | 10.95 | 17.20 | 24.07 | 44.60 | ▁▁▁▁▇ |
| Humedad_relativa | 0 | 1 | 39.48 | 51.22 | -200 | 34.05 | 48.55 | 61.88 | 88.73 | ▁▁▁▂▇ |
| Humedad_absoluta | 0 | 1 | -6.84 | 38.98 | -200 | 0.69 | 0.98 | 1.30 | 2.23 | ▁▁▁▁▇ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Fecha | 0 | 1 | 2004-03-10 | 2005-04-04 00:00:00 | 2004-09-21 00:00:00 | 391 |
| Hora | 0 | 1 | 1899-12-31 | 1899-12-31 23:00:00 | 1899-12-31 11:00:00 | 24 |
Con el resumen podemos observar que el conjunto de datos cuenta con
9357 observaciones de 15 variables.
R ha identificado dos tipos de datos en el
conjunto: numeric (15 variables) y
POSIXct (2 variables). Las variables de tipo
numeric corresponden a las mediciones de los sensores, mientras
que las dos variables de tipo POSIXct son la fecha y hora en
que se realizaron las mediciones.
Aunque parece que los datos están
completos, es importante recordar que los valores faltantes en el
conjunto fueron etiquetados con un valor de -200. Por
lo tanto, es necesario determinar la cantidad real de datos faltantes en
el conjunto tomando en cuenta esta etiqueta.
# Veamos un conteo de los datos faltantes para cada
# una de las variables.
sapply(datos_sucios, function(x) sum(x == -200))
## Fecha Hora
## 0 0
## Monoxido_carbono Oxido_estaño
## 1683 366
## Hidrocarburos_no_metanicos Benceno
## 8443 366
## Titanio Oxidos_nitrogeno
## 366 1639
## Oxido_tungsteno Dioxido_nitrogeno
## 366 1642
## Oxido_tungsteno_NO2 Oxido_indio
## 366 366
## Temperatura Humedad_relativa
## 366 366
## Humedad_absoluta
## 366
Ahora que hemos evaluado la cantidad real de datos faltantes en el conjunto de datos, podemos observar que la variable que contiene los datos de la concentración promedio de hidrocarburos no metánicos es la que tiene el mayor número de datos faltantes, con un total de 8443. Además, hay tres variables con más de 1600 datos faltantes y finalmente 9 variables con 366 datos faltantes, que corresponden a los datos obtenidos del sensor.
Para poder utilizar algunas de las funciones útiles de R destinadas específicamente a tratar con datos faltantes, es necesario cambiar los valores de -200 por NA. De esta manera, podremos aprovechar al máximo las herramientas disponibles en R para el análisis de datos.
# Reemplzamos los -200 por NA
datos_sucios <- replace(datos_sucios, datos_sucios == -200, NA)
# Usamos la libreria naniar para ver de forma grafica
# información de los datos faltantes
library(naniar)
## Warning: package 'naniar' was built under R version 4.3.2
# Porcentaje de NA
gg_miss_var(datos_sucios, show_pct = T)
# ¿En dónde se localizan los datos faltantes?:
vis_miss(datos_sucios)
De todas las variables disponibles en nuestro conjunto de datos nos centraremos en las siguientes variables, las cuales son fundamentales para la evaluación de la calidad del aire y son ampliamente monitoreadas en diferentes lugares del mundo:
Aunque las concentraciones de los compuestos químicos en la atmósfera son lo más relevante, consideramos la temperatura y la humedad relativa importantes para la evaluación de la calidad del aire, ya que pueden afectar la formación y dispersión de contaminantes en el aire.
Después de ver los datos faltantes y ver que la variable de la concentración promedio de hidrocarburos no metánicos es la que más datos faltantes tiene, la haremos a un lado para realizar la modelación con la mayor cantidad de datos.
# Seleccionamos las variables que usaremos para el
# análisis
library(dplyr)
calidad_Aire<- select(datos_sucios,
Monoxido_carbono,
Benceno,
Oxidos_nitrogeno,
Dioxido_nitrogeno,
Temperatura,
Humedad_relativa)
Finalmente, eliminemos todos los renglones donde exista al menos un dato faltante de nuestro conjunto de datos.
# Eliminamos renglones con CUALQUIER valor faltante
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.2
calidad_Aire_Limpio<- drop_na(calidad_Aire)
paste("La cantidad de observaciones que quedan después de eliminar los datos faltantes son: ", nrow(calidad_Aire_Limpio))
## [1] "La cantidad de observaciones que quedan después de eliminar los datos faltantes son: 6941"
Vemos que después de eliminar los renglones con al menos un dato faltante nos quedan nuestras 6 variables con 6941 observaciones; es decir, eliminamos aproximadamente el 26% de los datos.
Consideremos las variables que hemos seleccionado, es importante
estudiar sus distribuciones de probabilidad. Recordando un poco lo que
hicimos en la Fase 2 veamos si estas cumplen o no con
tener un perfil “normal”, veamos las distribuciones:
NOTA: La línea roja es la curva normal
teórica y la línea azul es la densidad empírica.
library(psych)
## Warning: package 'psych' was built under R version 4.3.2
library(paletteer) # Librería para las paletas de colores
## Warning: package 'paletteer' was built under R version 4.3.2
# Histogramas
multi.hist(x = calidad_Aire_Limpio, bcol=paletteer_d("RColorBrewer::BuPu"), dcol = c("blue", "red"), dlty = c("dotted", "solid"),lwd=c(2,3),
main = c("Monóxido de Carbono","Benceno",
"Óxidos de Nitrógeno","Dióxido de Nitrógeno",
"Temperatura","Humedad Relativa"), global = FALSE)
A primera impresión notamos que las variables correspondientes a las concentraciones de monóxido de carbono, benceno y óxidos de nitrógeno no parecen tener este perfil normal que estamos buscando, mientras que las variables que corresponden a la concentración de dióxido de nitrógeno, temperatura y humedad relativa parece que si cumplen con el perfil normal. Pero veamos unas graficas más:
attach(calidad_Aire_Limpio)
par(mfrow=c(2,3))
qqnorm(Monoxido_carbono, main="Normal Q-Q Plot (Monóxido de Carbono)",col="darkblue")
qqline(Monoxido_carbono, col="red", lwd=2)
qqnorm(Benceno, main="Normal Q-Q Plot (Benceno)",col="darkblue")
qqline(Benceno, col="red", lwd=2)
qqnorm(Oxidos_nitrogeno, main="Normal Q-Q Plot (Óxidos de Nitrógeno)",col="darkblue")
qqline(Oxidos_nitrogeno, col="red", lwd=2)
qqnorm(Dioxido_nitrogeno, main="Normal Q-Q Plot (Dióxido de Nitrógeno)",col="darkblue")
qqline(Dioxido_nitrogeno, col="red", lwd=2)
qqnorm(Temperatura, main="Normal Q-Q Plot (Temperatura)",col="darkblue")
qqline(Temperatura, col="red", lwd=2)
qqnorm(Humedad_relativa, main="Normal Q-Q Plot (Humedad Relativa)",col="darkblue")
qqline(Humedad_relativa, col="red", lwd=2)
Una forma de ver si existe normalidad consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta. Esta comparación corresponde a las graficas anteriores, con eso podemos concluir que nuestras sospechas eran ciertas, las observaciones de las variables correspondientes a las concentraciones de monóxido de carbono, benceno y óxidos de nitrógeno no caen sobre la recta lo que nos indica que no hay normalidad, para las otras tres variables podríamos decir que la “gran mayoría” de observaciones si caen sobre la recta. Incluso podríamos realizar una prueba mas rigurosa para comprobar la normalidad, pero por esta ocasión nos quedaremos y aceptaremos estas conclusiones. Entonces, tenemos tres variables problemáticas:
Concentración de monóxido de carbono
Concentración de benceno
Concentración de óxidos de nitrógeno
Profundicemos en detalle en estas tres variables
problemáticas y examinemos qué distribuciones de probabilidad
podrían modelarlas de manera más precisa.
A primera vista, las distribuciones de probabilidad de las tres variables problemáticas parecen exhibir un comportamiento similar (es curioso que tres variables sean descritas por una distribución normal y que las otras tres variables no tengan este comportamiento normal pero que sí sigan un comportamiento bastante parecido entre ellas). Si nos aventuramos aún más, podríamos afirmar que podemos modelar su comportamiento utilizando distribuciones Gamma. Sin embargo, sería prudente no aceptar esto como un hecho y analizar más detenidamente las variables:
Veamos primero la gráfica de Cullen &
Frey, esta gráfica se representa el coeficiente de
asimetría en el eje x y el coeficiente de curtosis en
el eje y. Cada punto en el gráfico representa una distribución de datos
específica. Si los datos siguen una distribución normal, los puntos se
agruparán cerca del (0,3) en el gráfico. Si los datos tienen una
asimetría o curtosis significativa, los puntos se alejarán de ese punto
y se desplazarán hacia diferentes direcciones.
La grafica de
Cullen & Frey nos ayuda a escoger alguno
de los modelos de probabilidad basados en la curtosis y la simetría.
# Creamos el grafico Cullen & Frey
library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 4.3.2
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: survival
descdist(data = Monoxido_carbono, graph = TRUE)
## summary statistics
## ------
## min: 0.1 max: 11.9
## median: 1.9
## mean: 2.182467
## estimated sd: 1.441158
## estimated skewness: 1.33915
## estimated kurtosis: 5.620881
Basado en el perfil del histograma y la gráfica de Cullen & Frey,
nuestras observaciones preliminares parecen confirmarse, y es plausible
que el modelo de probabilidad más apropiado sea una
Distribución Gamma, caracterizada por los
parámetros alfa y beta.
Veamos ahora el valor del
AIC (Akaike’s Information Criterion)
para diferentes distribuciones. El AIC es una medida
utilizada para la selección de modelos en el contexto del análisis de
datos y el ajuste de modelos estadísticos. Cuanto más bajo sea el valor
del AIC, mejor se considera que es el ajuste del
modelo. Esto significa que un modelo con un AIC más
bajo tiene un buen equilibrio entre ajuste y complejidad, y se prefiere
sobre modelos con AIC más alto.
library(univariateML)
## Warning: package 'univariateML' was built under R version 4.3.2
library(dplyr)
library(tibble)
## Warning: package 'tibble' was built under R version 4.3.2
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.
comparacion_aic <- AIC(
mlbetapr(Monoxido_carbono), # Beta prime distribution
mlexp(Monoxido_carbono), # Exponential distribution
mlinvgamma(Monoxido_carbono), # Inverse Gamma distribution
mlgamma(Monoxido_carbono), # Gamma distribution
mllnorm(Monoxido_carbono), # Log-normal distribution
mlrayleigh(Monoxido_carbono), # Rayleigh distribution
mlinvgauss(Monoxido_carbono), # Inverse Gaussian
mlweibull(Monoxido_carbono), # Weibull distribution
mlinvweibull(Monoxido_carbono) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>%
arrange(AIC)
## distribucion df AIC
## 1 mlgamma(Monoxido_carbono) 2 22420.80
## 2 mlweibull(Monoxido_carbono) 2 22602.19
## 3 mllnorm(Monoxido_carbono) 2 22810.90
## 4 mlbetapr(Monoxido_carbono) 2 22990.51
## 5 mlrayleigh(Monoxido_carbono) 1 23279.82
## 6 mlinvgauss(Monoxido_carbono) 2 23348.53
## 7 mlinvgamma(Monoxido_carbono) 2 24604.22
## 8 mlexp(Monoxido_carbono) 1 24718.29
## 9 mlinvweibull(Monoxido_carbono) 2 25344.70
Al analizar los valores del AIC de las diferentes
distribuciones, se observa que el valor más bajo corresponde a la
distribución Gamma, lo que indica que
efectivamente es la distribución más adecuada.
Hagamos el ajuste con
una distribución Gamma:
library(latex2exp)
## Warning: package 'latex2exp' was built under R version 4.3.2
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
library(paletteer)
#Realizamos el ajuste de la función de densidad Gamma
monoxido_fitGamma<- fitdistr(Monoxido_carbono, densfun="gamma")
# Graficamos el histograma junto con la curva de ajuste
hist(Monoxido_carbono,col=paletteer_d("ggsci::indigo_material"),
main=TeX("Monóxido de Carbono"),
xlab= TeX("$mg/m^3$"),
ylim = c(0, 0.4),
breaks=20,
prob=TRUE)
curve(dgamma(x, monoxido_fitGamma$estimate[1],
monoxido_fitGamma$estimate[2]),col="red",
lwd=2,add=T)
paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", monoxido_fitGamma$estimate[1], monoxido_fitGamma$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son: 2.34925510077524 1.07642238920912"
Para las otras dos variables problemáticas seguiremos el mismo procedimiento.
Veamos primero la gráfica de Cullen & Frey para la variable que contiene los datos de la concentración de Benceno:
# Creamos el grafico Cullen & Frey
library(fitdistrplus)
descdist(data = Benceno, graph = TRUE)
## summary statistics
## ------
## min: 0.1815254 max: 63.74148
## median: 8.788282
## mean: 10.55441
## estimated sd: 7.46517
## estimated skewness: 1.301695
## estimated kurtosis: 5.409269
Vemos que de acuerdo al perfil del histograma y la gráfica de
Cullen & Frey parece que nuestras sospechas se hacen
realidad, el modelo de probabilidad que parece más adecuado parece ser
nuevamente una Distribución Gamma de parámetros alfa y
beta.
Veamos ahora el valor del AIC
(Akaike’s Information Criterion) para diferentes
distribuciones.
library(univariateML)
library(dplyr)
library(tibble)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.
comparacion_aic <- AIC(
mlbetapr(Benceno), # Beta prime distribution
mlexp(Benceno), # Exponential distribution
mlinvgamma(Benceno), # Inverse Gamma distribution
mlgamma(Benceno), # Gamma distribution
mllnorm(Benceno), # Log-normal distribution
mlrayleigh(Benceno), # Rayleigh distribution
mlinvgauss(Benceno), # Inverse Gaussian
mlweibull(Benceno), # Weibull distribution
mlinvweibull(Benceno) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>%
arrange(AIC)
## distribucion df AIC
## 1 mlgamma(Benceno) 2 45042.35
## 2 mlweibull(Benceno) 2 45125.36
## 3 mllnorm(Benceno) 2 45583.01
## 4 mlinvgauss(Benceno) 2 46333.64
## 5 mlrayleigh(Benceno) 1 46411.88
## 6 mlexp(Benceno) 1 46597.54
## 7 mlbetapr(Benceno) 2 46711.92
## 8 mlinvgamma(Benceno) 2 47774.08
## 9 mlinvweibull(Benceno) 2 48181.68
Al analizar los valores del AIC de las diferentes
distribuciones, se observa que el valor más bajo corresponde a la
distribución Gamma, lo que indica que
efectivamente es la distribución más adecuada.
Hagamos el ajuste con
una distribución Gamma:
library(latex2exp)
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
library(paletteer)
#Realizamos el ajuste de la función de densidad Gamma
benceno_fitGamma<- fitdistr(Benceno, densfun="gamma")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
# Graficamos el histograma junto con la curva de ajuste
hist(Benceno,col=paletteer_d("ggsci::deep_orange_material"),
main=TeX("Benceno"),
xlab= TeX("$\\mu g/m^3$"),
ylim = c(0, 0.08),
breaks=20,
prob=TRUE)
curve(dgamma(x, benceno_fitGamma$estimate[1],
benceno_fitGamma$estimate[2]),col="red",
lwd=2,add=T)
paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", benceno_fitGamma$estimate[1], benceno_fitGamma$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son: 1.97409183844 0.187039581597446"
Estudiemos nuestra última variable problemática. Primero la gráfica de Cullen & Frey:
# Creamos el grafico Cullen & Frey
library(fitdistrplus)
descdist(data = Oxidos_nitrogeno, graph = TRUE)
## summary statistics
## ------
## min: 2 max: 1479
## median: 186
## mean: 250.6565
## estimated sd: 208.604
## estimated skewness: 1.643631
## estimated kurtosis: 6.155228
Para esta última variable, parece que la Distribución Gamma ya no es la más adecuada para ajustar los datos. Parece ser que una distribución beta podría ser la mejor opción. Sin embargo, antes de sacar conclusiones, es importante evaluar el valor del AIC (Akaike’s Information Criterion) para diferentes distribuciones.
library(univariateML)
library(dplyr)
library(tibble)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.
comparacion_aic <- AIC(
mlbetapr(Oxidos_nitrogeno), # Beta prime distribution
mlexp(Oxidos_nitrogeno), # Exponential distribution
mlinvgamma(Oxidos_nitrogeno), # Inverse Gamma distribution
mlgamma(Oxidos_nitrogeno), # Gamma distribution
mllnorm(Oxidos_nitrogeno), # Log-normal distribution
mlrayleigh(Oxidos_nitrogeno), # Rayleigh distribution
mlinvgauss(Oxidos_nitrogeno), # Inverse Gaussian
mlweibull(Oxidos_nitrogeno), # Weibull distribution
mlinvweibull(Oxidos_nitrogeno) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>%
arrange(AIC)
## distribucion df AIC
## 1 mllnorm(Oxidos_nitrogeno) 2 89673.24
## 2 mlgamma(Oxidos_nitrogeno) 2 89693.83
## 3 mlweibull(Oxidos_nitrogeno) 2 89884.58
## 4 mlinvgauss(Oxidos_nitrogeno) 2 90044.34
## 5 mlexp(Oxidos_nitrogeno) 1 90569.33
## 6 mlbetapr(Oxidos_nitrogeno) 2 91063.12
## 7 mlinvgamma(Oxidos_nitrogeno) 2 91137.15
## 8 mlinvweibull(Oxidos_nitrogeno) 2 91448.87
## 9 mlrayleigh(Oxidos_nitrogeno) 1 92914.62
Al analizar los valores del AIC de las diferentes
distribuciones, se observa que el valor más bajo corresponde a la
distribución Log-Normal, una distribución que
no teníamos en mente.
Hagamos el ajuste con una distribución
Log-Normal:
library(latex2exp)
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
library(paletteer)
#Realizamos el ajuste de la función de densidad Gamma
oxidos_fitLogNormal<- fitdistr(Oxidos_nitrogeno, densfun="lognormal")
# Graficamos el histograma junto con la curva de ajuste
hist(Oxidos_nitrogeno,col=paletteer_d("ggsci::teal_material"),
main=TeX("Óxidos de Nitrógeno"),
xlab= TeX("$ppb$"),
ylim = c(0, 0.0045),
breaks=20,
prob=TRUE)
curve(dlnorm(x, oxidos_fitLogNormal$estimate[1],
oxidos_fitLogNormal$estimate[2]),col="red",
lwd=2,add=T)
paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", oxidos_fitLogNormal$estimate[1], oxidos_fitLogNormal$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son: 5.18820287112231 0.862644837560372"
Entonces, podemos moodelar el comportamiento de la variable de
concentración de Óxidos de Nitrógeno con una distribución Log-Normal con
parametros:
meanlog = 5.188,
sdlog=0.862
Para transformar la(s) variable(s) para que tenga una distribución
normal usaremos la transformación de Box-Cox que es una
técnica estadística utilizada para estabilizar la varianza y mejorar la
normalidad de los datos. Esta técnica fue propuesta por los estadísticos
George Box y David Cox.
Para cada una de las variables problemáticas
realizaremos estas transformaciones y mostraremos a continuación los
histogramas y las gráficas Normal Q-Q para ver si la
transformación funcionó y las variables se normalizaron.
# Mónoxido de Carbono
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
library(paletteer)
library(latex2exp)
# Encontramos lambdas optimas
lambda_monoxido = BoxCox.lambda(Monoxido_carbono,method="loglik",
lower=-100,upper=100)
lambda_benceno = BoxCox.lambda(Benceno,method="loglik",
lower=-100,upper=100)
lambda_oxidos = BoxCox.lambda(Oxidos_nitrogeno,method="guerrero",
lower=-100,upper=100)
## Warning in optimize(guer.cv, c(lower, upper), x = x, nonseasonal.length =
## nonseasonal.length): NA/Inf replaced by maximum positive value
#Realizamos las transformaciones
trans.monoxido<- BoxCox(Monoxido_carbono, lambda_monoxido)
trans.benceno<- BoxCox(Benceno, lambda_benceno)
trans.oxidos<- BoxCox(Oxidos_nitrogeno, lambda_oxidos)
# Histogramas
# Monóxido de Carbono
par(mfrow=c(3,3))
hist(Monoxido_carbono,col=paletteer_d("ggsci::indigo_material"),
main=TeX("Monóxido de Carbono"),
breaks=20,
ylim = c(0, 0.4),
prob=TRUE) # variable original
hist(trans.monoxido, col=paletteer_dynamic("cartography::orange.pal", 20),
main=TeX("Transformación"),
breaks=20,
ylim = c(0, 0.5),
prob=TRUE) # variable transformada
# Graficas Q-Q
qqnorm(trans.monoxido, col="darkblue", cex = 0.7)
qqline(trans.monoxido, col="red", lwd=1)
# Benceno
hist(Benceno,col=paletteer_d("ggsci::deep_orange_material"),
main=TeX("Benceno"),
breaks=20,
ylim = c(0, 0.08),
prob=TRUE) # variable original
hist(trans.benceno, col=paletteer_dynamic("cartography::orange.pal", 20),
main=TeX("Transformación"),
breaks=20,
ylim = c(0, 0.3),
prob=TRUE) # variable transformada
# Graficas Q-Q
qqnorm(trans.benceno,col="darkblue", cex = 0.7)
qqline(trans.benceno, col="red", lwd=1)
# Óxidos de Nitrógeno
hist(Oxidos_nitrogeno,col=paletteer_d("ggsci::teal_material"),
main=TeX("Óxidos de Nitrógeno"),
breaks=20,
ylim = c(0, 0.0045),
prob=TRUE) # variable original
hist(trans.oxidos, col=paletteer_dynamic("cartography::orange.pal", 20),
main=TeX("Transformación"),
breaks=20,
ylim = c(0, 0.3),
prob=TRUE) # variable transformada
# Graficas Q-Q
qqnorm(trans.oxidos,col="darkblue", cex = 0.7)
qqline(trans.oxidos, col="red", lwd=1)
Podemos ver que después de la transformación realizada a las tres variables problemáticas, los perfiles de sus distribuciones (transformadas) y las gráficas Q-Q nos muestran aparente normalidad en todos los casos.
El monóxido de carbono (\(CO\)) es
un gas resultante de la combustión incompleta, que ocurre cuando existe
una baja concentración de oxígeno. Según la bibliografía,
aproximadamente el 86% de las emisiones de \(CO\) provienen del sector del transporte,
seguido por un 6% atribuido a la quema de combustibles en la industria,
y un 3% originado en procesos industriales. El 4% restante se genera a
través de quemas y otros procesos no identificados. Asimismo, el CO se
produce de manera natural mediante la oxidación del metano, un proceso
comúnmente asociado con la descomposición de materia orgánica.
El \(CO\) puede tener efectos
adversos en la salud, ya que compite con el oxígeno en la corriente
sanguínea, disminuyendo así la capacidad de transporte de oxígeno a los
diferentes órganos. Las personas sensibles, especialmente aquellas que
padecen problemas cardíacos, pueden experimentar una reducción en su
capacidad de oxigenación.
Dado lo anteriormente expuesto, resulta crucial llevar a cabo un
monitoreo de los niveles de CO. Por tanto, hemos decidido
seleccionar la concentración de monóxido de carbono (\(CO\)) como nuestra variable
respuesta.
Para realizar la división del conjunto de datos en dos grupos, se
propuso inicialmente basarla en la presencia de registros de
concentración de monóxido de carbono que excedieran los límites
establecidos por las organizaciones pertinentes. Por ejemplo, en México,
el límite para la concentración de monóxido de carbono como contaminante
atmosférico es de 12.595 miligramos por metro cúbico en un promedio
móvil de ocho horas al año, con el objetivo de proteger la salud de la
población susceptible. De manera similar, la Agencia de Protección
Ambiental de Estados Unidos (EPA, por sus siglas en inglés) ha
establecido un límite ambiental de 10 miligramos por metro cúbico para
el monóxido de carbono en el aire, promediado en un período de ocho
horas.
Sin embargo, surge un problema con esta división, ya que hay muy pocos registros que superen estas concentraciones. Si consideramos únicamente las observaciones que cumplen con un nivel superior a los 10 miligramos por metro cúbico, nos encontramos con tan solo 5 observaciones. Por lo tanto, se ha tomado la decisión de separar el conjunto de datos basándonos en la temperatura. Dado que las observaciones se obtuvieron en una ciudad italiana (no se especifica cuál), se han revisado los registros de temperatura promedio en la ciudad de Roma, Italia, correspondientes al año 2004, y se ha encontrado una temperatura media anual de 17.1°C. Al analizar los datos, se observa que el promedio de temperatura es de 17.7°C, lo que parece ser un criterio adecuado para la división. Además, resultará interesante comparar el promedio de la variable respuesta entre ambos grupos para ver si la concentración de monóxido de carbono aumenta o disminuye en función de la presencia de temperaturas por encima del promedio.
Hagamos entonces la separación:
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
datos_baja_temperatura<- calidad_Aire_Limpio[Temperatura <= 17.7, ]
datos_alta_temperatura<- calidad_Aire_Limpio[Temperatura > 17.7, ]
length(datos_baja_temperatura$Monoxido_carbono)
## [1] 3694
length(datos_alta_temperatura$Monoxido_carbono)
## [1] 3247
Vamos a evaluar si existe una diferencia significativa en el promedio
de la concentración de monóxido de carbono entre los dos grupos que
hemos creado. Para hacer esto, utilizaremos un intervalo de confianza
del 90%.
Inicialmente, observamos que la distribución del monóxido de carbono
no sigue una distribución normal. Sin embargo, para calcular el
intervalo de confianza del promedio, basta con tener mas de 25
observaciones para poder decir que los promedios se distribuyen
normales. Afortunadamente, contamos con 3694 observaciones en el grupo
datos_baja_temperatura y 3247 observaciones en el grupo
datos_alta_temperatura, por lo que se cumple este
requisito.
Dado que la varianza de la población es desconocida, utilizaremos un intervalo de confianza para saber si podemos asumir que son iguales.
# Veamos si es valido suponer varianzas iguales.
# Verifiquemos mediante un intervalo al 95%.
vartest.temp<- var.test(datos_alta_temperatura$Monoxido_carbono,
datos_baja_temperatura$Monoxido_carbono,
alternative = "two.sided",conf.level = 0.95)
vartest.temp$conf.int
## [1] 0.7410043 0.8467492
## attr(,"conf.level")
## [1] 0.95
Recordemos que con la prueba var.test lo que estamos viendo
es \(\sigma_1/\sigma_2\), si ambas
varianzas son iguales el intervalo de confianza de la prueba debe
contener al 1.
El intervalo de confianza calculado al 95% de
confianza es: \((0.741, 0.846)\), que
claramente no contiene al 1, por lo tanto no podemos asumir
varianzas iguales.
Finalmente podemos asumir que nuestras
observaciones son independientes. Encontremos entonces el intervalo de
confianza de la diferencia de promedios de la concentración de monóxido
de carbono entre los dos grupos que hemos creado.
# Realizamos la prueba
prueba_promedios<- t.test(datos_baja_temperatura$Monoxido_carbono,
datos_alta_temperatura$Monoxido_carbono,
alternative="two.sided",mu=0,paired=F,var.equal=F)
prueba_promedios$conf.int
## [1] -0.09108558 0.04382863
## attr(,"conf.level")
## [1] 0.95
Observamos que el valor cero se encuentra dentro de este
intervalo, además de ser un intervalo bastante estrecho con una longitud
de \(L = 0.134\). Por lo tanto,
llegamos a la conclusión de que no existen diferencias
significativas en el promedio de la concentración de monóxido de carbono
entre ambos grupos. Esto indica que la concentración promedio
de monóxido de carbono no experimenta cambios significativos en relación
a si las temperaturas registradas son superiores o inferiores a la
temperatura promedio.
Ahora procederemos a desarrollar un modelo que analice la posible
relación entre la concentración de monóxido de carbono y las
concentraciones promedio de benceno, \(NO_x\) y \(NO_2\), así como la influencia de la
temperatura y la humedad relativa en dichas concentraciones.
Antes
de proceder a la formulación de un modelo, es importante explorar más a
fondo las relaciones entre las distintas variables. A continuación,
analizaremos estas relaciones mediante representaciones gráficas:
library(GGally)
## Warning: package 'GGally' was built under R version 4.3.2
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.2
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
ggpairs(calidad_Aire_Limpio)
Veamos una gráfica más, una representación de la matriz de correlación:
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.2
## corrplot 0.92 loaded
corrplot(cor(calidad_Aire_Limpio))
De las gráficas presentadas podemos concluir que las variables que tienen una mayor relación lineal con la concentración de Monóxido de carboono son:
Concentración de Benceno (0.930).
Concentración de Óxidos de Nitrógeno (0.786).
Concentración de Dióxido de Nitrógeno (0.674).
Las variables predictoras que están más relacionadas entre sí son:
Óxidos de Nitrógeno - Benceno (0.718).
Dióxido de Nitrógeno - Benceno (0.604).
Dióxido de Nitrógeno - Óxidos de Nitrogéno (0.757).
Humedad Relativa - Temperatura (0.564).
Vamos a proceder a desarrollar un primer modelo que incluya todas las variables predictoras y evaluaremos su rendimiento inicial. Sin embargo, desde un primer vistazo, podemos anticipar que el modelo no será óptimo debido a la presencia de dependencia lineal entre varias variables predictoras. Además, debemos tener en cuenta que no estamos considerando en este modelo las variables que hemos transformado previamente y que las variables predictoras no siguen una distribución normal.
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 21):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores
model.monoxido1<- lm(Monoxido_carbono~Benceno+Oxidos_nitrogeno+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = calidad_Aire_Limpio)
summary(model.monoxido1)
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno +
## Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2598 -0.1851 0.0062 0.1922 3.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.836e-01 4.208e-02 4.364 1.3e-05 ***
## Benceno 1.566e-01 1.316e-03 119.015 < 2e-16 ***
## Oxidos_nitrogeno 8.093e-04 5.472e-05 14.788 < 2e-16 ***
## Dioxido_nitrogeno 2.465e-03 2.037e-04 12.099 < 2e-16 ***
## Temperatura -1.213e-02 9.522e-04 -12.736 < 2e-16 ***
## Humedad_relativa 1.585e-03 4.351e-04 3.643 0.000271 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4529 on 6935 degrees of freedom
## Multiple R-squared: 0.9013, Adjusted R-squared: 0.9012
## F-statistic: 1.267e+04 on 5 and 6935 DF, p-value: < 2.2e-16
Al parecer nuestra suposición inicial no era tan correcta, el primer modelo utilizando todas las variables predictoras no resulto tan malo, veamos algunas conclusiones:
El modelo con todas las variables introducidas como predictores tiene un R alto (0.90), es capaz de explicar el 90% de la variabilidad observada en la concentración de Monóxido de Carbono.
El valor-p del primer modelo es significativo (< 2.2e-16) por lo que se puede aceptar que el modelo no es por azar, al menos uno de los coeficientes parciales de regresión es distinto de 0.
Todas las variables son significativas (valor-p de todos los coeficientes del modelo).
Veamos un resumen del modelo:
library(gvlma)
valida_model <- gvlma(model.monoxido1)
summary(valida_model)
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno +
## Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.2598 -0.1851 0.0062 0.1922 3.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.836e-01 4.208e-02 4.364 1.3e-05 ***
## Benceno 1.566e-01 1.316e-03 119.015 < 2e-16 ***
## Oxidos_nitrogeno 8.093e-04 5.472e-05 14.788 < 2e-16 ***
## Dioxido_nitrogeno 2.465e-03 2.037e-04 12.099 < 2e-16 ***
## Temperatura -1.213e-02 9.522e-04 -12.736 < 2e-16 ***
## Humedad_relativa 1.585e-03 4.351e-04 3.643 0.000271 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4529 on 6935 degrees of freedom
## Multiple R-squared: 0.9013, Adjusted R-squared: 0.9012
## F-statistic: 1.267e+04 on 5 and 6935 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = model.monoxido1)
##
## Value p-value Decision
## Global Stat 2.193e+05 0.00000 Assumptions NOT satisfied!
## Skewness 2.013e+03 0.00000 Assumptions NOT satisfied!
## Kurtosis 2.168e+05 0.00000 Assumptions NOT satisfied!
## Link Function 4.058e+00 0.04397 Assumptions NOT satisfied!
## Heteroscedasticity 4.440e+02 0.00000 Assumptions NOT satisfied!
Como ya lo habíamos comentado, todas las variables son
significativas (valor-p de todos los coeficientes del
modelo) y deberían quedarse en el modelo.
Los errores estándar de
cada uno de los coeficientes son:
| Coeficiente | Error Estándar |
|---|---|
| (Intercept) | 4.208e-02 |
| Benceno | 1.316e-03 |
| Oxidos_nitrogeno | 5.472e-05 |
| Dioxido_nitrogeno | 2.037e-04 |
| Temperatura | 9.522e-04 |
| Humedad_relativa | 4.351e-04 |
Y el error estándar residual del modelo es:
0.4529.
Diagnóstico mediante las gráficas que proporciona el programa R
# Graficas de diagnostico
plot(model.monoxido1)
Veamoos algunoos comentarios respecto a las gráficas.
Gráfica residuos vs valores ajustados: Los residuos parecen estar concentrandose en la primera parte de la gráfica, mostrando alguna especie de patrón. Las observaciones 5817, 4258 y 4259 presentan diferencias mayores, por lo que quizá estas observaciones “tengan algún problema”.
Normalidad residuos: Los residuos parecen no comportarse con normalidad, y vemos los mayores problemas en las “colas”.
Gráfica localización-escala: Vemos que la varianza de los residuos ocomienza a crecer, formando una especie de cono. La línea roja va en aumento y no es hoorizoontal.
Gráfica residuos vs leverage: El residuo 4259 está fuera de la distancia de Cook.
Lo que parecía un buen modelo por el alto valor de \(R^2\) vemos que tiene problemas muy importantes, no está cumpliendo ninguno de los supuestos para poder hacer la regresión lineal, veamos si podemos mejorarlo un poco.
step(object = model.monoxido1, direction = "both", trace = 1)
## Start: AIC=-10989.16
## Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + Dioxido_nitrogeno +
## Temperatura + Humedad_relativa
##
## Df Sum of Sq RSS AIC
## <none> 1422.6 -10989.2
## - Humedad_relativa 1 2.72 1425.3 -10977.9
## - Dioxido_nitrogeno 1 30.03 1452.6 -10846.2
## - Temperatura 1 33.28 1455.9 -10830.7
## - Oxidos_nitrogeno 1 44.86 1467.5 -10775.7
## - Benceno 1 2905.65 4328.3 -3268.1
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno +
## Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
##
## Coefficients:
## (Intercept) Benceno Oxidos_nitrogeno Dioxido_nitrogeno
## 0.1836346 0.1566329 0.0008093 0.0024648
## Temperatura Humedad_relativa
## -0.0121272 0.0015850
# Agregamos los datos al conjunto
calidad_Aire_Limpio$transformacion_benceno <- trans.benceno
calidad_Aire_Limpio$transformacion_oxidos <- trans.oxidos
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 23):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores transformados
model.monoxido2<- lm(Monoxido_carbono~transformacion_benceno+transformacion_oxidos+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = calidad_Aire_Limpio)
summary(model.monoxido2)
##
## Call:
## lm(formula = Monoxido_carbono ~ transformacion_benceno + transformacion_oxidos +
## Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1199 -0.3180 -0.0677 0.2559 5.4020
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5265823 0.0666031 -7.906 3.06e-15 ***
## transformacion_benceno 0.8537410 0.0103556 82.443 < 2e-16 ***
## transformacion_oxidos 0.0394882 0.0091440 4.318 1.59e-05 ***
## Dioxido_nitrogeno 0.0014394 0.0003096 4.649 3.40e-06 ***
## Temperatura -0.0263575 0.0014257 -18.487 < 2e-16 ***
## Humedad_relativa 0.0014952 0.0005825 2.567 0.0103 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.597 on 6935 degrees of freedom
## Multiple R-squared: 0.8285, Adjusted R-squared: 0.8284
## F-statistic: 6701 on 5 and 6935 DF, p-value: < 2.2e-16
# Graficas de diagnostico
plot(model.monoxido2)
Ya habíamos dicho que uno de los supuestos que debe satisfacer el modelo de regresión múltiple es el de no multicolinealidad, para ver si existe usaremos: Tolerancia y Factor de Inflación de la Varianza (VIF), si:
VIF > 10 es causa de preocupación.
VIF es sustancialmente mayor que 1 entonces la regresión puede verse perjudicada.
Tolerancia = 1/VIF debajo de 0.1 indica un problema serio.
Tolerancia debajo de 0.2 indica un problema potencial.
# Criterio VIF
library(car)
## Warning: package 'car' was built under R version 4.3.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.2
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
vif(model.monoxido1)
## Benceno Oxidos_nitrogeno Dioxido_nitrogeno Temperatura
## 3.265584 4.408959 3.163125 2.399888
## Humedad_relativa
## 1.946206
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura, transformacion_benceno,
## transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 18):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 20):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 26):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores transformados
model.monoxido3<- lm(Monoxido_carbono~Benceno+Dioxido_nitrogeno+Temperatura, data = calidad_Aire_Limpio)
summary(model.monoxido3)
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Dioxido_nitrogeno +
## Temperatura, data = calidad_Aire_Limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7656 -0.1933 0.0026 0.1911 3.1530
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3489067 0.0213727 16.32 <2e-16 ***
## Benceno 0.1709388 0.0010206 167.49 <2e-16 ***
## Dioxido_nitrogeno 0.0034282 0.0001614 21.25 <2e-16 ***
## Temperatura -0.0203306 0.0007031 -28.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4619 on 6937 degrees of freedom
## Multiple R-squared: 0.8973, Adjusted R-squared: 0.8973
## F-statistic: 2.021e+04 on 3 and 6937 DF, p-value: < 2.2e-16
library(gvlma)
valida_model <- gvlma(model.monoxido3)
summary(valida_model)
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Dioxido_nitrogeno +
## Temperatura, data = calidad_Aire_Limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7656 -0.1933 0.0026 0.1911 3.1530
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3489067 0.0213727 16.32 <2e-16 ***
## Benceno 0.1709388 0.0010206 167.49 <2e-16 ***
## Dioxido_nitrogeno 0.0034282 0.0001614 21.25 <2e-16 ***
## Temperatura -0.0203306 0.0007031 -28.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4619 on 6937 degrees of freedom
## Multiple R-squared: 0.8973, Adjusted R-squared: 0.8973
## F-statistic: 2.021e+04 on 3 and 6937 DF, p-value: < 2.2e-16
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = model.monoxido3)
##
## Value p-value Decision
## Global Stat 2.719e+05 0.000 Assumptions NOT satisfied!
## Skewness 2.010e+03 0.000 Assumptions NOT satisfied!
## Kurtosis 2.693e+05 0.000 Assumptions NOT satisfied!
## Link Function 1.874e+00 0.171 Assumptions acceptable.
## Heteroscedasticity 5.628e+02 0.000 Assumptions NOT satisfied!
# Graficas de diagnostico
plot(model.monoxido2)
boxplot.matrix(as.matrix(calidad_Aire_Limpio),use.cols = T)
library(olsrr)
## Warning: package 'olsrr' was built under R version 4.3.2
##
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
##
## cement
## The following object is masked from 'package:datasets':
##
## rivers
ols_plot_cooksd_chart(model.monoxido1)
# DFBETA mide la diferencia en la estimacion de cada parametro con y sin el punto influyente.
# Hay un valor de DFBETA para cada observacion, es decir, hay n observaciones y k parametros, por lo que hay n*k DFBETAs.
# En general, valores grandes de DFBETAS indican observaciones que son influyentes en la estimacion de cierto parametro.
ols_plot_dfbetas(model.monoxido1)
ols_plot_resid_stud(model.monoxido1)
library(car)
outlierTest(model.monoxido1)
## rstudent unadjusted p-value Bonferroni p
## 4259 -18.844505 2.7147e-77 1.8843e-73
## 4258 -14.065135 2.5194e-44 1.7487e-40
## 5817 -7.434707 1.1740e-13 8.1486e-10
## 4918 6.648333 3.1897e-11 2.2140e-07
## 4204 6.577777 5.1231e-11 3.5559e-07
## 5418 -6.274431 3.7204e-10 2.5824e-06
## 4260 -6.111312 1.0417e-09 7.2304e-06
## 5417 -6.103445 1.0940e-09 7.5936e-06
## 5406 -5.940388 2.9811e-09 2.0692e-05
## 5407 -5.816380 6.2812e-09 4.3598e-05
registros_a_eliminar <- c(4259, 4258, 5817, 4918, 4204, 5418, 4260, 5417, 5406, 5407)
datos_sin_registros <- calidad_Aire_Limpio[-registros_a_eliminar, ]
model.monoxido4<- lm(Monoxido_carbono~Benceno+Oxidos_nitrogeno+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = datos_sin_registros)
summary(model.monoxido4)
##
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno +
## Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = datos_sin_registros)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.58037 -0.18460 0.00603 0.19348 2.56970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.314e-01 3.949e-02 5.861 4.82e-09 ***
## Benceno 1.619e-01 1.250e-03 129.578 < 2e-16 ***
## Oxidos_nitrogeno 7.070e-04 5.143e-05 13.747 < 2e-16 ***
## Dioxido_nitrogeno 2.345e-03 1.910e-04 12.278 < 2e-16 ***
## Temperatura -1.468e-02 8.972e-04 -16.366 < 2e-16 ***
## Humedad_relativa 1.270e-03 4.079e-04 3.114 0.00186 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.424 on 6925 degrees of freedom
## Multiple R-squared: 0.9131, Adjusted R-squared: 0.913
## F-statistic: 1.455e+04 on 5 and 6925 DF, p-value: < 2.2e-16
plot(model.monoxido4)
Después de eliminar las observaciones influyentes, hemos observado
una mejora en el valor de \(R^2\) en
nuestro modelo, lo cual puede parecer positivo a simple vista. Sin
embargo, al examinar detenidamente las demás métricas y las últimas
gráficas, es evidente que el modelo no es sólido. A pesar de que
tradicionalmente utilizamos el valor de \(R^2\) como una medida para determinar la
calidad del modelo, en este caso nos encontramos con problemas que
incumplen los supuestos necesarios para realizar la regresión. Ninguno
de estos supuestos se cumple, lo cual es motivo de preocupación y nos
lleva a cuestionar si este modelo puede ser utilizado de manera
apropiada.
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura, transformacion_benceno,
## transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura, transformacion_benceno,
## transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 18):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 19):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 20):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 22):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 28):
##
## Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
## Oxidos_nitrogeno, Temperatura
nuevos_datos <- data.frame(
var1 = c(max(Benceno)+5,
max(Oxidos_nitrogeno)+20,
max(Dioxido_nitrogeno)+20,
max(Temperatura)+3,
max(Humedad_relativa)+2), # Valores para var1
var2 = c(max(Benceno)+10,
max(Oxidos_nitrogeno)+40,
max(Dioxido_nitrogeno)+40,
max(Temperatura)+5,
max(Humedad_relativa)+4), # Valores para var2
var3 = c(max(Benceno)+15,
max(Oxidos_nitrogeno)+50,
max(Dioxido_nitrogeno)+50,
max(Temperatura)+7,
max(Humedad_relativa)+6), # Valores para var3
stringsAsFactors = FALSE
)
predicciones <- predict(model.monoxido1, nuevos_datos, interval="prediction")
## Warning: 'newdata' had 5 rows but variables found have 6941 rows
head(predicciones)
## fit lwr upr
## 1 2.370100 1.4820364 3.258164
## 2 1.879968 0.9918811 2.768055
## 3 1.921226 1.0331577 2.809294
## 4 2.030767 1.1426972 2.918838
## 5 1.555741 0.6676693 2.443813
## 6 1.193149 0.3050997 2.081198